% This code goes through various levels of disaggregation, based on the
% actual inputs for price stickiness, sector size and I/O linkages
% and produces reduced sized input according to sector size

% by Raphael Schoenle, August 10, 2016 and April 2019

clearvars -except wd
cd(fullfile(wd, 'Figure4_RHS')); 
%clc; 
%clear all;
%cd '/Users/raphaelschoenle/Dropbox/Networks/code/JME_post/Figure4_RHS'

%% initial data
% Data
load FPA.mat
%Parameters values
params.rho_a = 0.9;         % persistence agg tech shocks
params.rho_mu = 0.9;         % persistence monetary shock
params.rho_k = 0.9;         % persistence idiosync tech shock
params.beta = 0.9975;      % discount factor
params.theta = 6;          % elast of subst within sectors
params.delta = 0.5;

Ck = csvread("2002CkDetailed.csv");
Znum = csvread("2002RevShareDetailed.csv");
n = length(Ck); % Number of Industries
Z = Znum(1:n,1:n);
Wc = Ck/sum(Ck); % Share of consumption; the weights omega

%number of sectors
n_sectors = 341;
n_sectors1 = 58;

%number of periods in IRF
n_periods = 30;
n_periods_idio = 30;

%number of cases
num_case = 18;

% consumption weights
Omega_C_het = Wc;
Omega_C_vec_pareto = Omega_C_het./sum(Omega_C_het);

% sectoral Calvo
load FPA.mat
FPA_het = FPA;
alpha_vec_pareto = FPA_het;

%sectoral IO weights 
W = Z./(ones(n,1)*sum(Z)); %Normalize the share of inputs for each sector, omega
psi = params.delta*(params.theta - 1)/params.theta;
psi_C = (1-params.delta*(params.theta - 1)/params.theta);
n_k = inv(eye(n)-psi.*W)*(psi_C*Wc);
Omega_het = W;
Omega_IO_matrix_pareto = W;

% I/O matrix in terms of values (total value of shipments provided by a supplier)
total_inputs = xlsread('intermediate_use.xlsx');

cd(fullfile(wd, 'Figure4_RHS')); 
%cd '/Users/raphaelschoenle/Dropbox/Networks/code/JME_post/Figure4_RHS'

% produce aggregation up to 341-n_reductions sectors only
n_reductions = 340;

% set weights to 0 when aggregating frequencies
Omega_C_het_W = Omega_C_het;
Omega_C_het_W(Omega_C_het_W<0)=0;
Omega_C_het_W=Omega_C_het_W./sum(Omega_C_het_W);

for jj=1:n_reductions
    
    Omega_C_het_sorted = sort(Omega_C_het);
    dalpha = log(Omega_C_het_sorted(2:end)./Omega_C_het_sorted(1:end-1));
    dalpha_min = min(dalpha);

    alpha_merge_pos = find(dalpha_min==dalpha);

    % first, if there is identical size of sectors: 
    % merge first such two together (notation in terms of alpha)
    if dalpha_min == 0
        merge_pos = find(Omega_C_het_sorted(alpha_merge_pos(1))==Omega_C_het);
        
        % new data
        Omega_C_het_new = [Omega_C_het(1:merge_pos(1)-1); ...
            Omega_C_het(merge_pos(1))+Omega_C_het(merge_pos(2));... 
            Omega_C_het(merge_pos(1)+1:merge_pos(2)-1);...
            Omega_C_het(merge_pos(2)+1:end)];
        Omega_C_het_W_new = [Omega_C_het_W(1:merge_pos(1)-1); ...
            Omega_C_het_W(merge_pos(1))+Omega_C_het_W(merge_pos(2));... 
            Omega_C_het_W(merge_pos(1)+1:merge_pos(2)-1);...
            Omega_C_het_W(merge_pos(2)+1:end)];
        
        % aggregate up total input usage
        total_inputs_new = [total_inputs(1:merge_pos(1)-1); ...
            total_inputs(merge_pos(1))+total_inputs(merge_pos(2));... 
            total_inputs(merge_pos(1)+1:merge_pos(2)-1);...
            total_inputs(merge_pos(2)+1:end)];
        
        if ( Omega_C_het_W(merge_pos(1))...
            +Omega_C_het_W(merge_pos(2))) == 0
        FPA_het_new = [FPA_het(1:merge_pos(1)-1);...
            0.5*FPA_het(merge_pos(1))+...
            0.5*FPA_het(merge_pos(2));...
        FPA_het(merge_pos(1)+1:merge_pos(2)-1);...
        FPA_het(merge_pos(2)+1:end)];
        else
            FPA_het_new = [FPA_het(1:merge_pos(1)-1);...
            Omega_C_het_W(merge_pos(1))/( Omega_C_het_W(merge_pos(1))...
            +Omega_C_het_W(merge_pos(2)))*FPA_het(merge_pos(1))+...
            Omega_C_het_W(merge_pos(2))/( Omega_C_het_W(merge_pos(1))...
            +Omega_C_het_W(merge_pos(2)))*FPA_het(merge_pos(2));...
        FPA_het(merge_pos(1)+1:merge_pos(2)-1);...
        FPA_het(merge_pos(2)+1:end)];
        end;
        
        Omega_het = Omega_het';
        % add columns
            Omega_het_new = [Omega_het(:,1:merge_pos(1)-1) ...
            Omega_het(:,merge_pos(1))+Omega_het(:,merge_pos(2)) ...
            Omega_het(:,merge_pos(1)+1:merge_pos(2)-1)...
            Omega_het(:,merge_pos(2)+1:end)];
        
        % add rows
        Omega_het_new = [Omega_het_new(1:merge_pos(1)-1,:); ...
            total_inputs(merge_pos(1))/( total_inputs(merge_pos(1))+ total_inputs(merge_pos(2)))*Omega_het_new(merge_pos(1),:)+...
            total_inputs(merge_pos(2))/( total_inputs(merge_pos(1))+ total_inputs(merge_pos(2)))*Omega_het_new(merge_pos(2),:);...
            Omega_het_new(merge_pos(1)+1:merge_pos(2)-1,:);...
            Omega_het_new(merge_pos(2)+1:end,:)];

        Omega_het_new = Omega_het_new';

        % renormalize I/O matrix so that sum across columns, given row, is
        % equal to 1
        Omega_het_new = Omega_het_new./(ones(size(Omega_het_new,1),1)*sum(Omega_het_new)); %Normalize the share of inputs for each sector, omega
    else
        merge_pos1 = min(find(Omega_C_het_sorted(alpha_merge_pos(1))==Omega_C_het), find(Omega_C_het_sorted(alpha_merge_pos(1)+1)==Omega_C_het));
        merge_pos2 = max(find(Omega_C_het_sorted(alpha_merge_pos(1))==Omega_C_het), find(Omega_C_het_sorted(alpha_merge_pos(1)+1)==Omega_C_het));
         
        if ( Omega_C_het_W(merge_pos1)...
            +Omega_C_het_W(merge_pos2)) == 0
        FPA_het_new = [FPA_het(1:merge_pos1-1);...
            0.5*FPA_het(merge_pos1)+...
            0.5*FPA_het(merge_pos2);...
        FPA_het(merge_pos1+1:merge_pos2-1);...
        FPA_het(merge_pos2+1:end)];
        else
            FPA_het_new = [FPA_het(1:merge_pos1-1);...
            Omega_C_het_W(merge_pos1)/( Omega_C_het_W(merge_pos1)...
            +Omega_C_het_W(merge_pos2))*FPA_het(merge_pos1)+...
            Omega_C_het_W(merge_pos2)/( Omega_C_het_W(merge_pos1)...
            +Omega_C_het_W(merge_pos2))*FPA_het(merge_pos2);...
        FPA_het(merge_pos1+1:merge_pos2-1);...
        FPA_het(merge_pos2+1:end)];
        end
                  
        Omega_C_het_new = [Omega_C_het(1:merge_pos1-1); ...
            Omega_C_het(merge_pos1)+Omega_C_het(merge_pos2);... 
            Omega_C_het(merge_pos1+1:merge_pos2-1);...
            Omega_C_het(merge_pos2+1:end)];
        
        Omega_C_het_W_new = [Omega_C_het_W(1:merge_pos1-1); ...
            Omega_C_het_W(merge_pos1)+Omega_C_het_W(merge_pos2);... 
            Omega_C_het_W(merge_pos1+1:merge_pos2-1);...
            Omega_C_het_W(merge_pos2+1:end)];

        % aggregate up total input usage
        total_inputs_new = [total_inputs(1:merge_pos1-1); ...
            total_inputs(merge_pos1)+total_inputs(merge_pos2);... 
            total_inputs(merge_pos1+1:merge_pos2-1);...
            total_inputs(merge_pos2+1:end)];

       Omega_het = Omega_het';
       
       % add columns
            Omega_het_new = [Omega_het(:,1:merge_pos1-1) ...
            Omega_het(:,merge_pos1)+Omega_het(:,merge_pos2) ...
            Omega_het(:,merge_pos1+1:merge_pos2-1)...
            Omega_het(:,merge_pos2+1:end)];
        
        % add rows
        Omega_het_new = [Omega_het_new(1:merge_pos1-1,:); ...
            total_inputs(merge_pos1)/( total_inputs(merge_pos1)+ total_inputs(merge_pos2))*Omega_het_new(merge_pos1,:)+...
            total_inputs(merge_pos2)/( total_inputs(merge_pos1)+ total_inputs(merge_pos2))*Omega_het_new(merge_pos2,:);...
            Omega_het_new(merge_pos1+1:merge_pos2-1,:);...
            Omega_het_new(merge_pos2+1:end,:)];
        
        Omega_het_new = Omega_het_new';
        
        % renormalize I/O matrix so that sum across columns, given row, is
        % equal to 1
        Omega_het_new = Omega_het_new./(ones(size(Omega_het_new,1),1)*sum(Omega_het_new)); %Normalize the share of inputs for each sector, omega

    end

    FPA_het=FPA_het_new;
    total_inputs = total_inputs_new;
    Omega_C_het=Omega_C_het_new;
    Omega_C_het_W=Omega_C_het_W_new;
    Omega_het=Omega_het_new;
    filename = ['agg_data_' num2str(n_sectors-jj) '_C.mat'];
    save(filename, 'FPA_het', 'Omega_C_het' , 'Omega_het', 'Omega_C_het_W');
end

% the following part makes a plot of the distributions of alphas, and
% computes the variance across aggregation levels

clearvars -except wd
%clear all
n_reductions = 340;
n_red=[ 1 141 241 291 316 331 333];

cd(fullfile(wd, 'Figure4_RHS')); 
%cd '/Users/raphaelschoenle/Dropbox/Networks/code/JME_post/Figure4_RHS'

filename = ['agg_data_' num2str(341-n_red(1)) '_C.mat']
load(filename)

    
    % variance
	clear FPA FPA_het Omega_C Omega
    n_reductions = 340
    for i_k=1:n_reductions
        filename = ['agg_data_' num2str(341-i_k) '_C.mat'];
        load(filename);
        vfpa(i_k)=var(FPA_het);
    end
    
    figure()
    counts = (341-1):-1:(341-n_reductions);
    line=plot( counts, vfpa, 'linewidth',2)
    xlabel('\textbf{Number of sectors}','Interpreter','latex')
    ylabel('\textbf{Variance of frequency of price changes}','Interpreter','latex')

	NameArray = {'LineStyle', 'Color','Marker'}
     ValueArray = {'-','b','none'}
     set(line ,NameArray, ValueArray(1,:))

     filename = ['variance_C_aggregation.pdf'];
    saveas(gca,filename);
    
    % distributions of modified first order outdegrees
	for ll = 1:n_reductions
        filename = ['agg_data_' num2str(341-ll) '_C.mat'];
        load(filename);
        fo_outdegrees(1:length(Omega_C_het),ll) = (Omega_C_het.*(1-FPA_het))'*Omega_het';
    end
            
        
